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We further develop the theory of quantum finite-size effects in metallic nanoparticles, which was 
originally formulated by Ifache, Ricard and Flytzanis [J. Opt. Soc. Am. B 3, 1647 (1986)] and 
(in a somewhat corrected form) by Rautian [Sov. Phys. JETP 85, 451 (1997)]. These references 
consider a metal nanoparticle as a degenerate Fermi gas of conduction electrons in an infinitely-high 
spherical potential well. This model (referred to as the HRFR model below) yields mathematical 
expressions for the linear and the third-order nonlinear polarizabilities of a nanoparticle in terms 
of infinite nested series. These series have not been evaluated numerically so far and, in the case 
of nonlinear polarizability, they can not be evaluated with the use of conventional computers due 
to the high computational complexity involved. Rautian has derived a set of remarkable analytical 
approximations to the series but direct numerical verification of Rautian's approximate formulas 
remained a formidable challenge. In this work, we derive an expression for the third-order nonlinear 
polarizability, which is exact within the HRFR model but amenable to numerical implementation. 
We then evaluate the expressions obtained by us numerically for both linear and nonlinear polariz- 
abilities. We investigate the limits of applicability of Rautian's approximations and find that they 
are surprizingly accurate in a wide range of physical parameters. We also discuss the limits of 
small frequencies (comparable or below the Drude relaxation constant) and of large particle sizes 
(the bulk limit) and show that these limits are problematic for the HRFR model, irrespectively of 
any additional approximations used. Finally, we compare the HRFR model to the purely classical 
theory of nonlinear polarization of metal nanoparticles developed by us earlier [Phys. Rev. Lett. 
100, 47402 (2008)]. 



I. INTRODUCTION 

This paper is dedicated to the memory of Sergey Gle- 
bovich Rautian (1928-2009) who was a teacher to some 
of us and inspiration to all. 

Metal nanoparticles have received an extraordinary 
amount of attention recently because of their ability to 
greatly enhance local fields. The enhancement is at- 
tributed to the excitation of surface plasmons and it 
has a variety of applications in photovoltaicai, sensing^ 
and surface-enhanced Raman scattering^-^^. Currently, 
nanopartices of very small sizes, up to a few nanometers, 
are customarily used in experiments. The theoretical de- 
scription of the optical properties of such nanoparticles 
is most frequently based on the macroscopic electrody- 
namics. At least, this is typical in the field of plasmonics. 
However, macroscopic electrodynamics can not capture 
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certain effects of finite size. Hache, Ricard and Flytza- 
nis^ and Rautian (in a somewhat modified form)^- have 
developed an elaborate theory of quantum finite-size ef- 
fects in metal nanoparticles. In Refs. [^0, a nanoparticle 
was modeled as degenerate Fermi gas confined in an infi- 
nite potential well of spherical shape (below, the HRFR 
model). Despite being fairly simple, the HRFR model re- 
sults in very complicated formulas, which can not be eval- 
uated numerically even with the aid of modern comput- 
ers. For example, the expression for the third-order non- 
linear polarizability involves a twelve-fold nested sum- 
mation. Rautian has reduced the number of summations 
from twelve to eight by performing summation over the 
magnetic sublevels analytically; he then obtained a num- 
ber of remarkable approximations to the resulting eight- 
fold summation^. However, these approximations have 
never been verified directly due to the overwhelming nu- 
merical complexity involved. 

In this contribution, we develop the analytical the- 
ory of Rautian a step further by reducing the number 
of nested summations involved from eight to five with- 
out making any additional approximations. This turns 
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out to be sufficient to render the formulas amenable to 
direct numerical implementation. We then compare the 
resuhs of numerical evaluation of the five-fold series de- 
rived by us to the resuhs, which fohow from Rautian's ap- 
proximate formulas, and discuss various physical limits, 
including the limits of low frequency and large particle 
size. 

In Sec. In] we review Rautian's theory. Here we use 
somewhat simplified notation and, in particular, avoid 
the use of irreducible spherical tensors and 6j-symbols. 
In Sec, mil we develop the theory further by utilizing the 
orbital selection rules and reduce the nested summation 
involved in the definition of the third-order nonlinear po- 
larizability from eight-fold to five-fold. In Sec. IIV| we 
describe a simple method to relate the internal and ap- 
plied fields, which is to the first order consistent with the 
approach proposed in Ref. |^, but is more mathematically 
rigorous. In Sec. [Vj the results of numerical computa- 
tions are reported. A summary of obtained results and a 
discussion are contained in Sec. I VII 



II. RAUTIAN'S THEORY 

We start by reviewing Rautian's theory of quantum 
finite-size effects in conducting nanoparticles^. The 
physical system under consideration is a gas of ,jY 
non-interacting electrons placed inside a spherical, in- 
finitely deep potential well of radius a and subjected to 
harmonically-oscillating, spatially-uniform electric field 



Ei(t) = Ai exp(— iwt) -I- c.c. 



(1) 



Note that is the electric field inside the nanoparticle. 
It will be related to the external (applied) field Eg in 
Sec. HV] below. 

Since the nanoparticle is assumed to be electrically 
small (that is, a ^ A = l-KcjuS), the electron-field in- 
teraction can be described in the dipole approximation 
by the time-dependent operator 



V{t) = d ■ E(i) , 



(2) 



where d is the dipole moment operator. Under the addi- 
tional assumption that Ei(i) is linearly polarized with a 
purely real amplitude A^ 



AiZ, we can write 



V{t) = Gexp(-«wt) -|-( 



where 



G 



d ■ A, 



(3) 



(4) 



Rautian made use of the interaction representation in 
which the wave function is expanded in the basis of 
the unperturbed Hamiltonian eigenstates. The single- 
electron unperturbed states are 



<^nZm(r) 



1 
Znl 



(5) 



where ji{x) are the spherical Bessel functions of the 
first kind and order I; is the n-th positive root 
(n — 1,2,.. .) of the equation ji{x) — 0, Yim{i^) are the 
spherical functions (viewed here as functions of the polar 
and azimuthal angles of the unit vector f — r/r) and 



Ji+l(Cn 



(6) 



are normalization factors. The energy eigenstates are la- 
beled by the main quantum number n, the orbital num- 
ber I and the magnetic number m. The unperturbed 
energy levels are given by the formula 



where 



Enl — Eo^„i , 



En = 



(7) 



(8) 



and me is the electron mass. 

In what follows, we use the composite indices /i, v, 
77 and C to label the eigenstates. Each composite index 
corresponds to the triplet of quantum numbers {n,l,m). 
By convention, if /i = (n, I, m), then fi' = (n', Z', m'). The 
matrix elements of the z-projection of the dipole moment 
operator are given by 



z • d„„/ = eaAaa' , 



(9) 



where 



Smm'Rnl (bimSl-lJ' + bi+i^rnSl+l,l') , (10) 



Sw are the Kronecker delta-symbols, 1,1' > and 

'iS.nlS.n'l' 



km = 



nn I 

^nl — 



4/2-1 ' (.2 -|2 

\'inl 'in' I') 



(11) 



Note that the diagonal elements of A are all equal to zero, 
as is the case for any system with a center of symmetry. 
Finally, the matrix elements of the operator G are given 
by 



(12) 



The density matrix of the system, p, can be written in 
the form p^^(t) = p^^(t) exp(ia;^^i), where = {E^^ - 
E^)/h are the transition frequencies and Pfi^{t) is the so- 
called slow-varying amplitude, which obeys the following 
master equation^: 



d_ 



V 



(13) 



Here A^^ are the equilibrium state populations and F^,^ 
are phenomenological relaxation constants. Following 
Rautian, we assume that 



= Fi(5py + F2(l - Sf^iy) 



(14) 
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Eq. (|14p is the least complex assumption on F^i^, which 
still distinguishes the relaxation rates for the diagonal 
and off-diagonal elements of the density matrix. 

It can be seen that, for the case of zero external field, 
P/ii/ = S^uN^. The Fermi statistics is introduced at this 
point by writing 

2 

" exp[(£;^ - EF)l{kBT)] + 1 ' ^^^^ 

where Ep is the Fermi energy, ks is Boltzmann's con- 
stant, T is the temperature, and the factor of two in the 
numerator accounts for the electron spin. Conservation 
of particles reads = jY . When ^ > 1, the 

well-known analytical formula for the Fermi's energy. 



2\2/3 



(16) 



holds with a good accuracy. Here ^ is the characteristic 
atomic scale, defined by the relation 



where 



1^ = Vt/^ 



n = 47raV3 



(17) 



(18) 



is the nanoparticle volume. Thus, is the specific vol- 
ume per conduction electron. We note that £ is, gener- 
ally, different from the lattice constant h. Many met- 
als of interest in plasmonics have an fee lattice struc- 
ture with four conduction electrons per unit cell. In this 
case h — A^^'^i. For example, in silver, £ « 0.26nm, 
h « 0.41nm and Ep ~ 5.51eV. At room temperature 
(T = 300K), ksT « 0.026eV, so that T = is a good 
approximation. In this case, iV^ = 2 if E"^ < Ep and 
Nfj^ — otherwise. Most numerical results shown below 
have been obtained in this limit. However, to illustrate 
the effects of finite temperature, we have also performed 
some computations at T = 300K. Finally, the Fermi 
velocity is given by the equation 



Vp 



2Ep 



(19) 



In silver, vp « 1.2 • lO^cm/sec and, correspondingly, 
c/vp « 250. Electron velocities in excited states are 
expected to be no larger than a few times Fermi veloc- 
ity, still much smaller than c. This justifies the use of 
non-rclativistic quantum mechanics. 

The solution to (11311 has the form of a Fourier series: 



(20) 



Piiuit)= ^ exp(-isa;i) 



s— — oo 



The expansion coefficients piil obey the system of equa- 
tions: 



- G ( + 



(21) 



where 



iV 



(22) 



are Lorentzian spectral factors. 

The optical response of the nanoparticle is deter- 
mined by the quantum-mechanical expectation of its to- 
tal dipole moment, which is given by 



{d(t)) = ea^ A^,,p^i.(i) . 



(23) 



Upon substitution of ([20]) into ([SS]) . we obtain the ex- 
pansion of {d{t)) into temporal Fourier harmonics. We 
now consider the optical response at the fundamental fre- 
quency Lo, which describes degenerate nonlinear phenom- 
ena such as the four- wave mixing. Denoting the compo- 
nent of {d{t)), which oscillates at the frequency lu, by 
(duiit)), we can write 



where 



D 



/it/ 



(24) 



(25) 



The coefficients and the amplitude D in ([24|) can be 
expanded in powers of Ai. Namely, we can write 



D = nxiA,)A, 



(26) 



where 



XiA) =X1+ X3iA/A,t? + X5(A.Mat)" . . . (27) 

Here we have introduced the characteristic atomic field 



A^t = e/£^ 



(28) 



and have used the assumption that Ai is real-valued; in 
the more general case, the expansion contains the terms 
of the form Xsl^i/^atP, etc. Note that the definition of 
X3 in (I26l) . (|27p is somewhat unconventional. The nonlin- 
ear susceptibility x^^\ defined in standard expositions 
of the subject ^S, has the dimensionality of the inverse 
square of the electric field, that is, of cm^/statvolt^ = 
cm'^/erg in the Gaussian system of units. Here we find 
it more expedient to define dimensionless coefficients xij 
X3j X5j etc., and expand the dipole moment amplitude 
D in powers of the dimensionless variable Ai/A^t- 

The first two coefficients in the expansion (l?71) . xi and 
X3, have been computed by Rautian explicitly and are 
given by the following series: 



Xi 



X3 



(ea)4A2 



(29a) 



) 

at 



J2 S^IJA^cAc^A^.A,^ , (29b) 
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where 



A(o) 



1/ \ r/u 



A 



A«iV,c(AW+A(-)) {AfJ 



A(o) 

fJ-V 



(30) 



Here A^^j, = — N^. Note that all quantities inside 
the summation symbols are dimensionless and so are the 
factors in front of the summation symbols. 

The expression (|29bp involves a staggering 12-fold sum- 
mation (recall that each composite index fi, v, rj and ( 
consists of three integer indices) . Rautian used the math- 
ematical formalism of irreducible spherical tensors and 
6j-symbols to perform summation over magnetic sub- 
levels analytically and to reduce the expression to an 8- 
fold summation. However, this approach does not make 
use of the orbital selection rules, which are explicit in 
([TU| . In Sec. IHI[ we will use the orbital selection rules 
to analytically reduce (|29bp to a 5-fold summation. The 
resultant formula is amenable to direct numerical imple- 
mentation, as will be illustrated in Sec. |Vl 

Having performed the summation over the magnetic 
sublevels, Rautian has evaluated the resulting series by 
exploiting the following two approximations: 

1. Adopt the two- level approximation^. In this ap- 
proximation, only the terms with fJ. = rj and v = C 
are retained in the right-hand side of (j29b[) . 

2. Assume that there are two dominant contributions 
to the series (pUj) . The off-resonant (Drudean) con- 
tribution is obtained by keeping only the terms 
with Wpi, <C w in the Lorentzian factors A|fj(w). 
The resonant contribution is obtained by keeping 
only the terms with w^i/ ~ lu. Each contribution is 
then evaluated separately by replacing summation 
with integration. 

Using the same approximations, we have reproduced 
Rautian's analytical results. For xi, we obtained 



1 



Xi 



p 



Fi 



vp/a 



where 



' 47re2 



(31) 



(32) 



is the plasma frequency and Fi , gi are dimensionless real- 
valued functions, which weakly depend on the parame- 
ters of the problem and are of the order of unity. More 
specifically, Fi is very close to unity for all reasonable 
particle sizes and approaches unity asymptotically when 
a — >• oo (we have verified this numerically). In what fol- 
lows, we assume that _Fi = 1. The function gi depends 



most profoundly on the ratio k = fno/Ep. We can write, 
approximately. 



1 



51 



(33) 



An analytical expression for this integral and a plot are 
given in the Appendix. 

Equation ((3T|) is equivalent to combining equations 
3.16 and 3.23 of Ref. 0. On physical grounds, one 
can argue that these expressions are applicable only if 
T2/L1J ^ 1. Indeed, in the classical Drudc model, we 
have 



Drudc 
Al 



An uj(lo + i"f) 



(34) 



where 7 is a relaxation constant. We expect the Drude 
model to be accurate in the limit a — >■ c», when the sec- 
ond term in the square brackets in pip vanishes. Thus, 
(1311) has an incorrect low-frequency asymptote. We ar- 
gue that the asymptote is incorrect because the HRFR 
model disregards the Hartree interaction potential. This 
will be discussed in more detail in Sec. IIVI below. At 
this point, we assume that T2/UJ ^ 1 and expand (PT|) 
in T2/UJ, neglect the correction to the real part of the 
resulting expression, and obtain: 



Xi 



At: \ uj J 



1 



.2r2 + givp/a 



(35) 



This expression corresponds to equation 3.28 of Ref. 0. 
Note that neglecting the correction to the real part but 
retaining the correction to the imaginary part in the 
above equation is mathematically justified because the 
terms 2r2 and vp/a can be of the same order of magni- 
tude, as we will see below. 

Comparing (j35p to a similar expansion of (1341) , we con- 
clude that the size-dependent relaxation constant 7 is 
given by 



7 « 7oo + 51 



Vp 



7o 



2r2 



(36) 



where 700 is the relaxation constant in bulk. It can be 
seen that the ratio vp/a plays the role of the collision 
frequency. The analytical result (j36]) is widely known and 
frequently used; it will be confirmed by direct numerical 
evaluation of (|29ap below. 

For the third-order nonlinear susceptibility X3, we ob- 
tain, with the same accuracy as above. 



X3 



2 fay 
TiOttS \1) 



T2 a 



A. 



F. 



F. 



gs- 



(37) 



Here a = /He ~ 1/137 is the fine structure constant, 
Ap — 2'Kc/LOp is the wavelength at the plasma frequency 
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(ss 138nm in silver) and F^, is another set of dimen- 
sionless real- valued functions of the order of unity. For re- 
alistic parameters, the function F3 varies only slightly^iS 
between 0.30 and 0.33; we have taken F3 = 0.33 in the 
numerical computations of Sec. |Vl The function 173 can 
be approximated by the following integral: 

53«- r x^/'ix + nf/'dx, ^=^. (38) 

1^ Jl-K 

The approximate formula (j38|) applies only for hw < Ep. 
However, we are interested in the spectral region a; < cjp. 
In silver, fujjp sa 8.98eV and Hujp/Ep w 1.63. This leaves 
us with the spectral range Ep/h < uj < uip in which ([38)) 
is not applicable. The integral (|38p can be evaluated 
analytically; the resulting expression and plot are given 
in the Appendix. 

Expression (l37|) contains several dimensionless param- 
eters. For a silver nanoparticle of the radius a = lOnm, 

IOttS \eJ \£ J 

The ratio r2/ri is more puzzling. While r2 can be re- 
lated to the Drude relaxation constant through ([36|) . Fi 
does not enter the analytical approximations (pit . ([55)1 
or the exact expression (|29ap . Therefore, Fi can not be 
directly related to any measurement of the linear opti- 
cal response. It was previously suggested^ that, based 
on the available experimental studies of non-equilibrium 
electron kinetics in silver—"—, F2/F1 ^ 10. This ratio 
will be employed below. 

Another interesting question is the dependence of the 
results on the particle radius, a. It follows from the an- 
alytical approximation (|35|) that xi approaches a well- 
defined "bulk" limit when a — >■ 00. The characteristic 
length scale is vp/'joo (~ 44nm in silver). Of course, 
direct numerical evaluation of xi according to (|29ap is 
expected to reveal some dependence of xi on a, which is 
not contained in the analytical approximation (1351) . and 
this fact will be demonstrated below in Sec. IV Bl How- 
ever, it will also be demonstrated that (|35|) becomes very 
accurate in the spectral range of interest when a > 5nm. 
Thus, the HRFR model yields a result for xii which is 
consistent with the macroscopic limit. 

The situation is dramatically different in the case of 
the nonlinear susceptibility xa- It follows from p7l) that 
-^3 °~^°°> 0{a^). Therefore, there is no "bulk" limit for 
X3. This is an unexpected result. While some studies 
suggest that a positive correlation between xs s^nd a in a 
limited range of a is consistent with experimental mea- 
surements^, we can not expect this correlation to hold 
for arbitrarily large values of a, as this would, essentially, 
entail an infinite value of xa in bulk. Such a prediction 
appears to be unphysical. Of course, Rautian's theory 
is not expected to apply to arbitrarily large values of a 
because the interaction potential ([2]) is written in the 
dipole approximation and, moreover, it assumes that the 
electric field inside the nanoparticle is potential, that is. 



V X E = is a good approximation. Still, the absence 
of a "bulk" limit for X3 is troublesome. We, therefore, 
wish to understand whether the quadratic dependence of 
X3 on a is a property of the HRFR model itself or an ar- 
tifact of the additional approximations made in deriving 
the analytical expression ([37| . More specifically, we can 
state the following two hypotheses: 

1. The quadratic dependence of xs on a is an artifact 
of the approximations made in deriving the ana- 
lytical expression ([57]) from (j29b[) [these approxi- 
mations are listed explicitly between Eqs. ([30)) and 
([3T|) ]. In this case, we can expect that direct evalu- 
ation of (j29b[) will not exhibit the quadratic growth. 

2. The quadratic dependence of X3 on a is a prop- 
erty of the HRFR model itself. In particular, the 
absence of a "bulk" limit for xa can be caused by 
the following reasons: (i) The HRFR model ne- 
glects the retardation effects in large particles, (ii) 
The HRFR model does not account for the Hartree 
interaction potential. (In reality, interaction of 
the conduction electrons with the induced charge 
density may be important, especially, for comput- 
ing nonlinear corrections.) (iii) The HRFR model 
makes use of a phenomenological boundary condi- 
tion at the nanoparticle surface. 

Verification of these hypotheses was previously hindered 
by the computational complexity of the problem. In what 
follows, we render Rautian's theory amenable to direct 
numerical validation. Then we show that the analytical 
approximation ([37]) is surprisingly good. Therefore, the 
second hypothesis must be correct. 



III. RAUTIAN'S THEORY FURTHER 
DEVELOPED 

It is possible to simplify ([29bp without adopting any 
approximations. To this end, we deviate from Rautian's 
approach of using irreducible spherical tensors and 6j- 
symbols. Instead, we directly substitute the expressions 
([lD,([Tn|),([IIl) into ()29bl) . We use the selection rules in 
([T0| and the following results: 

' 4 _ IjiP + l) 
= - i5(4P-l) ' (^^^^ 

S, ^ ± >iA.,, - (39b) 

m— — l 

to evaluate summations over all magnetic quantum num- 
bers and over all orbital quantum numbers but one. This 
leaves us with a five-fold summation over four main quan- 
tum numbers and one orbital quantum number. After 
some rearrangements, we arrive at the following expres- 
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variables: 



X3 



J2 {ZiPi + SiQi) 



1=1 



where 



P, 



E 



B 



713,/ — 1,714,/ 



b: 



,712,/ — l,n4,/ 

Til ,/,n3 ,/ — 1 



,ri3,(— 1 



ni ,/,n2 ,/+l 
na,/ — 1,714,/ 



n2,/+l,n4,/ 
ni ,/,n3 — 1 



ni,n2,n3,'i4 



I r,n4,,l,n3,l — l 



B. 



ns — l,ni ,/ 
n4,/,n2,/H-l 



pn2,i+l pn4,' pn3,'-l 



(40) 



(41a) 



(41b) 



This expression is exact within the HRFR model. The 
two-level approximation corresponds to keeping only the 
first term in the brackets in (1401) and, further, keeping 
only the terms with n2 = and ni = in (j41ap . 



D = nA^t (xix + X3x\xf + X5x\x\ 



D = 
where 

X 



f^^at (aiy + a3y\y\ 



a5y\y\ 



D 



An 



y 



= y 



[aiy + a3y\y\ 



(42a) 
(42b) 



(43) 



Here we have accounted for the fact that there can be a 
phase shift between the internal and the external fields; 
therefore, Ai and A^. can not be real-valued simultane- 
ously. In the theory presented above, we assume that 
Ai is real- valued, and this can always be guaranteed by 
appropriately choosing the time origin. In this case, Ag 
is expected to be complex. 

The coefficients Xk in (|42ap can be found from 
Rautian's theory; our task is to find the coefficients 
in (I42bp given the constraint (|43l) . To this end, we sub- 
stitute (|43p into (|42ap and obtain a series in the variable 
y. We then require that the coefficients in this series and 
in (j42bp coincide. This yields an infinite set of equations 
for ak , the first two of which read 



IV. RELATING THE INTERNAL AND THE 
APPLIED FIELDS 

In the HRFR model, electrons move in a given, 
spatially- uniform internal field ([T]). In practice, one is 
interested in the optical response of the nanoparticle to 
the external (applied) field. We denote the amplitude of 
the external field by Ag = A^z. The two fields differ 
because of a charge density induced in the nanoparti- 
cle. The interaction of the conduction electrons with the 
induced charge density is described by the Hartree po- 
tential. However, rigorous introduction of the Hartree 
interaction into the HRFR model is problematic. Doing 
so would require the mathematical apparatus of density- 
functional theory. We can, however, apply here the clas- 
sical concept of the depolarizing field, although this ap- 
proach is less fundamental. 

In the macroscopic theory, a sphere (either dielectric 
or conducting) , when placed in a spatially- uniform, qua- 
sistatic electric field of frequency w and amplitude Ae, is 
polarized and acquires a dipole moment of an amplitude 
D. The electric field inside the sphere is also spatially- 
uniform and has the amplitude A;. The induced charge 
accumulates at the sphere surface in a layer whose width 
is neglected. Under these conditions, A^ = Ag — D/a'^. 
Note that a linear dependence between D and Ag is not 
assumed here. The form of the depolarizing field, — D/a"^, 
follows only from the assumption of spatial uniformity of 
the internal field and from the usual boundary conditions 
at the sphere surface. Then the Hartree interaction can 
be taken into account as follows. 

Let us introduce the dimensionless variables x = 
Ai/A^t and y = Ag/A^t. Then we can expand D in both 



Xi 



X3 



4:77 

1 a-i 

3 



ai 



in 
~3 



■ai 



An 
"3 



■ai 



47r 



"3 



(44a) 
(44b) 



It is convenient to introduce the linear field enhancement 
factor /i according to 



/i 



1 



1 + (47r/3)xi 



(1 



(45) 



where ei = 1 + 47rxi is the linear dielectric permittivity. 
Then the solutions to (l44t have the form 



ai=/iXi , a3 = /i'l/i|'x3 • (46) 

The factor /, which relates the external and internal field 
amplitudes according to Ai = fA^, is then found from 



X 4n , 

y 3 



a3\y\ 



(47) 



Using pS)) . we find that, to first order in ///at = lyP = 

l^eMatP, 



/ = ./l-^/?l/l|^X3f 



(48) 



Here we have introduced the intensity of the incident 
beam, / = (c/27r)|v4ep, and the "atomic" intensity /at = 

(c/2^)|^atp. 

Note that our approach to finding the field enhance- 
ment factor / is somewhat different from that adopted in 
Ref. 8, where the expansions (|42al) has been truncated at 
the third order and the truncated expression was assumed 
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to be exact. The results obtained in the two approaches 
coincide to first order in ///at- In Ref;^, higher order 
corrections to this result have also been obtained. In our 
approach, these corrections depend on the higher-order 
coefficients X5, xii ^tc, which have not been computed 
by Rautian. 

We finally note that the phenomenological accounting 
for the Hartree interaction described in this section, while 
is necessary for comparison with the experiment, does 
not remove the two main difficulties of the HRFR model. 
Specifically, it does not fix the low-frequency limit for xi 
and does not affect the oc a? dependence of X3- Regard- 
ing the low- frequency limit, we note that limt^_j.o / = 
and the internal field in the nanoparticle tends to zero 
in this limit. The induced macroscopic charge is local- 
ized at the sphere surface where the electric field jumps 
abruptly. In a more accurate microscopic picture, the 
width of this surface layer is finite and the electric field 
changes smoothly over the width of this layer. Unfortu- 
nately, the classical concept of depolarizing field can not 
capture surface phenomena of this kind. 
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FIG. 1: (color online) Energy eigenstates, which enter the 
summation according to (|40p .(|41 [) . for a = lOnm and /max = 
= 120. The horizontal blue line shows the Fermi energy. 
In this example, the total number of states below the Fermi 
surface is .yV ~ 2.4 • 10^ (counting all degeneracies) and the 
total number of states shown in the figure is 2nniax/max(/max4- 
2) ^3.5- 10^ 



V. NUMERICAL RESULTS 
A. Convergence 

We have computed the Bessel function zeros using 
the method of bisection and achieved a numerical dis- 
crepancy of the equation ji{^ni) = of less than 10~^^ 
for all values of indices. Since the function ji{x) is ap- 
proximately linear near its roots, we believe that we have 
computed with sufficiently high precision. 

The summation over I in pO)) was truncated so that 
I < 'max and the quadruple summation in (|4T|) was trun- 
cated so that ni, 712, na, 714 < ^^max- A typical set of en- 
ergy levels used in the summation is shown in Fig. [1] for 
the case a = lOnm, /max = '^max = 120. Here the en- 
ergy levels (normalized to Eq) are shown by dots and the 
horizontal axis corresponds to the orbital number, /. Re- 
ferring to Fig. [U we note that /max has been chosen so 
that all states with / > /max are above the Fermi surface. 
Since the electron transitions occur between two states 
with / and /' such that /' = / ± 1, the factors A^^^ for 
any transition involving the states with / > /max are zero 
(or exponentially small at finite temperatures). It can be 
seen that convergence with / is very fast - contribution of 
the terms in (PCT|) with / > /max is either zero (at T = 0) 
or exponentially small. 

The choice of 7imax is a more subtle matter. Since there 
are no selection rules on n, transitions can occur between 
two states (one below and one above the Fermi surface) 
with very different values of n and, correspondingly, very 
different energies. However, transitions with energy gaps 
much larger than fko are suppressed by the Lorentzian 
factors (j22p . In most numerical examples, we have cho- 
sen rimax so as to account for, at least, all transitions 
with the energy gaps of AE < 3hid. Many (but not all) 



transitions with larger energy gaps were also accounted 
for. This approach yields a result with seven significant 
figures. However, it results in too many terms in the 
summation when huj ~ Ep and a > lOnm. For these 
values of parameters, we have used a smaller rimax so as 
to account for, at least, all transitions with AE < hto. 
We estimate that the relative error incurred by this trun- 
cation is < 10%. 



B. Linear response 

We begin by considering the linear susceptibility xi- In 
computations, we use the commonly accepted parameters 
for silver, hujp = 8.98eV (Ap — 2TTc/ujp = 138nm) and 
7oo/wp = 0.002. Here the relaxation constant which 
enters (I29ap . is determined from = 2r2 (see ([55)) ). 
The frequencies used satisfy the condition 700 <C 1. 
More specifically, the ratio uj/ujp varies in the range 
0.02 < uj/ujp < 1. We do not consider the frequencies 
above ujp because silver exhibits strong interband absorp- 
tion in that spectral range. Except when noted other- 
wise, all computations have been carried out at T = 0. 

Fig. [2] displays the quantity — Rexi computed numer- 
ically by direct evaluation of (I29ap and by the Drude 
formula (j34p with the size-corrected relaxation constant 
7 ([551) . The factor gi in ([55)) has been computed using 
the analytical formula (|A1[) . At sufficiently high frequen- 
cies, the Drude model predicts that — 47rRexi ~ i^p/^)'^ 
and this behavior is reproduced for all radiuses consid- 
ered with good precision. However, at smaller frequen- 
cies, there are differences between the analytical approx- 
imation and the numerical results. These differences are 
especially apparent for a = 2nm. In this case, the opti- 
cal response of the sphere is, effectively, dielectric rather 



than metallic for uj < O.OBwp. A similar behavior has 
been observed at a = 4nm (data not shown). The emer- 
gence of a dielectric response in metal nanoparticles of 
sufficiently small size at sufficiently low frequencies has 
been overlooked in the past. It occurs due to discreteness 
of the energy states. Consider a particle with a = 2nm 
at zero temperature. In this case, the lowest-energy elec- 
tronic transition, which is allowed by Fermi statistics 
(that is, a transition with N^i, 7^ 0), occurs between the 
states {n = 1,1 = 18) and (n' 1,/' = 19). The corre- 
sponding transition frequency is uJmin ~ 0.056ajp. It can 
be seen from Fig.[2^a) that the particle becomes dielectric 

for UJ < Wmin- 

We now turn to consideration of the relaxation phe- 
nomena. To this end, we plot in Fig. [3] the quantity 



as a function of frequency. We note that Z, as defined in 
(j49| . is positive for all passive materials and, in the Drudc 
model, Z = 7/0;^; here 7 is size-corrected. It can be seen 
that the analytical formula (I36|) captures the relaxation 
phenomena in the nanoparticle surprisingly well. How- 
ever, as in Fig.[2Ia), the analytical approximation breaks 
down when a — 2nm and uj < O.OGwp. A similar break- 
down was observed for a — 4nm (data not shown) . For all 
other values of parameters, the numerically computed Z 
is reasonably close to the size-corrected value of 7/wp and 
exhibits the same overall behavior. The small systematic 
error at higher frequencies is, most likely, caused by the 
approximation (1331) for gi. It was, in fact, mentioned by 
Rautian that (j33p is hardly accurate when huj ^ Ep. 

The fine structure visible in Fig. [31[a,b) is due to dis- 
creteness of electron states. The allowed transition fre- 
quencies can be "grouped" , which results in the ap- 
pearance of somewhat broader peaks, clearly seen in 
Fig.|3Ub,c). While spectral signatures of discrete states in 
metal nanoparticles have been observed experimentally 
(including the effect of "grouping" the positions of in- 
dividual spectral peaks should not be invested with too 
much significance. In any realistic system, these peaks 
will be smoothed out by particle polydispersity, varia- 
tions in shape, and by nonradiative relaxation and energy 
transfer to the surrounding medium. 

The finite-size correction (|36)) to the Drude relaxation 
constant is widely known and used. However, the deriva- 
tions of have been, so far, either heuristic or relied 
on poorly controlled approximations. In Fig. [31 we have 
provided, to the best of our knowledge, the first direct, 
first- principle numerical verification of (|36p and of its lim- 
its of applicability. 

C. Nonlinear response 

We next turn to the nonlinear susceptibility X3- The 
same parameters for silver as before will be used. In 
addition, the calculations require the relaxation constant 
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FIG. 2: (color online) The quantity — Re^i as a function of 
frequency for particles of different radius a, as labeled. Cen- 
tered symbols correspond to direct numerical evaluation of 
(|29a|l and continuous curves show the Drude formula (|34p in 
which the size-corrected relaxation constant 7 (|36|) has been 
used. 



Fi. As was mentioned above, the experimental value of 
Fi can not be inferred by observing the linear optical 
response. It was previously suggested^ that F2/ri « 10. 
This value will be used below. 

In Fig. m we plot the absolute value of X3 as a function 
of the particle radius a for u = O.lup and for the Frohlich 
frequency oj = uip/y/S ~ 0.58wp, and compare the results 
of direct numerical evaluation of (j40|) to the analytical 
approximation ([37]). In the case uj — O.lwp, the analytical 
approximation is very accurate for a > 8nm and gives 



9 



1 F 



0.1 



0.01 ^ 



0.001 



Z Numerical 
(a) a = 2nm /, , 



0.1 F 



0.01 - 



0.001 



0.1 F 



■ z 



0.01 - 



0.001 




0.1 



Numerical 
(c) a — 16nm 



0.1 



Lj/ujp 



1 1 1 1 1 1 1 1 


Numerical : 

(b) a = 8nm ^/^p ■ 

Jco/LOp ■ 








0.1 t^/^^p 1 




FIG. 3: (color online) The quantity Z defined in (|49|l as a 
function of frequency for different particle radiuses, as labeled. 
Results of direct numerical computation are compared to the 
size-corrected Drude relaxation constant 7 (given by (|36|l ') and 
to its bulk value 700 (obtained in the limit a — > 00). 



in this case and results in a series of sharp maxima and 
minima of 1x3(^^)1- is shown in the inset of Fig. UJ^a), 
the function |x3((j)| is discontinuous. These discontinu- 
ities are artifacts of the zero temperature approximation. 
The introduction of a finite temperature (T — 300K) re- 
moves the discontinuities (see the inset) but does not 
eliminate the fine structure of the curve. Note, however, 
that the computations have been carried out with a very 
fine step in a, which is, arguably, unphysical: the param- 
eter a in a real nanosphere can change only in quantized 
steps of the order of the lattice constant h (« 0.41nm 
in silver). Moreover, the fine structure of X3 is unlikely 
to be observable experimentally due to the unavoidable 
effects of particle polydispersity. Therefore, the general 
trend given by the analytical approximation p7p can be 
a more realistic estimate of X3 for a < 4nm. 

Next consider the large-a behavior. For a > 8nni, the 
analytical and the "exact" formulas yield results, which 
are scarcely distinguishable. In particular, the quadratic 
growth of X3 with a has been confirmed up to a = 64nm 
in the case w = O.lwp - the largest radius for which nu- 
merical evaluation of (|40p is still feasible. This confirms 
Hypothesis 2 stated above, namely, that the quadratic 
growth of X3(«) is a property of the HRFR model itself 
rather than of the additional approximations, which were 
made to derive the analytical results. 

In Fig. O we study the dependence of |x3| on the fre- 
quency uj for fixed values of a. It can be seen that the 
accuracy of Rautian's approximation improves for larger 
particles and higher frequencies. At a = lOnm, the ap- 
proximation is nearly perfect in the full spectral range 
considered. 



D. Magnitude of the nonlinear effect and 
comparison with the classical theory of electron 
confinement 

In the previous subsection, we have plotted the coef- 
ficient X3j which appears in the expansion (|27p. The 
dimensionless parameter of this expansion, Ai/A^t, con- 
tains the amplitude of the internal electric field, Ai. How- 
ever, it is the amplitude of the external (applied) field, 
Ae, which can be directly controlled in an experiment. 
The incident beam intensity is given by / = (c/27r)|^eP- 
We can use the results of Sec. llVlto write 



the correct overall trend for a < 4nm. A systematic 
discrepancy of unknown origin between the approximate 
and the numerical results is observed for 4nm < a < 8nni. 
In the case uj = ujp/y/3, the analytical approximation 
gives the correct trend in the whole range of a considered. 
Note that, in the case u! = O.lcjp, the absolute value of 
X3 is dominated by Rex3 and Imxs for large and small 
values of a, respectively. When uj = Wp/VS, the real part 
of X3 is dominating for all values of a used in the figure. 

Consider first particles with a < 4nm. As expected, 
the discreteness of energy levels plays an important role 



D = nA, [ai+az{I/h 



(50) 



where ai and are related to xi and X3 by (j46p and 
-^at = (c/27r)A^j is the characteristic "atomic" intensity. 
The quantity /at can be expressed in terms of the funda- 
mental physical constants and the material-specific pa- 
rameter In the case of silver, /at ~ 2.3 • 10^'*W/cm^. 
Obviously, intensities of such magnitude are not achiev- 
able in any experiment. However, the magnitude of the 
nonlinear correction can be amplified by the two impor- 
tant effectsiS,: the effect of synchronism (constructive in- 
terference), which is not considered here, and by the ef- 
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FIG. 4: (color online) Absolute value of the nonlinear sus- 
ceptibility, 1x3 1, computed by direct evaluation of (I40|l and 
by analytical approximation (|37|) as a function of the particle 
radius for cj = O.ltJp (a) and for ui = uip /^/3 ^ 0.58tjp (b). 
Logarithmic scale is used on both axes. The inset in panel 
(a) shows a zoom of the plot for 3.4nm < a < 3.7nm. In the 
inset, the results of evaluating (|40[) at T = and at T = 3007^ 
are shown. 



feet of local field enhancement, which will be taken into 
account by using the expressions derived in Sec. IIVI 

We will also compare the expression (fSOl) to the re- 
sults obtained from the purely classical argumentsi^. In 
Ref. [l^, we have argued the surface charge in a polarized 
metal nanoparticle can not be confined to an infinitely 
thin layer. When the width of this layer is not negligible 
(compared to the particle radius) , a nonlinear correction 
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FIG. 5: (color online) Absolute value of the nonlinear suscep- 
tibility, 1x3 1, computed by direct evaluation of (|40p at zero 
temperature and by analytical approximation (|37|) as a func- 
tion of the frequency for a — 4nm (a) and for a = lOnm 
(b). 



to the particle polarizability is obtained. After some re- 
arrangement of the formulas derived in this reference, we 
can express the amplitude D, defined analogously to (|24|) . 
as 



D = nA, 



/3i+/33 07^+-. 



where 



/3i 



3 

An ujp/3 — up- — i^uj 



33 = /3i|/3i| 

TT a 



(51) 

(52a) 
(52b) 



We note that ai > /3i. That is, the linear po- 

larizabilities of both theories are the same in the region 
of parameters where the theories are applicable. The 
classical theory, however, does not contain the low fre- 
quency anomaly in the linear polarizability. On the other 
hand, relaxation is introduced in Ref . [H through the phe- 
nomenological parameter 7 whose dependence on a can 
not be deduced theoretically. Below, we use the result of 
the quantum theory, namely, Eq. p6l) for the relaxation 
constant 7 (|52ap . 
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It can be seen that the classical and quantum expres- 
sions for D are quite different. The first non- vanishing 
nonhnear correction in (j52b|) is of the order of ^///at 
but contains an additional small parameter £/a. Thus, 
the nonlinear correction depends differently on the in- 
cident intensity, frequency and the particle radius a in 
the two theories. Additionally, Rautian's theory con- 
tains the parameter T1/T2, which does not enter into 
the classical theory. These factors complicate a direct 
comparison of the two results. We will, therefore, fo- 
cus on the trends for one particular value of the incident 
power, / — lOkW/cm^. One should bear in mind that 
the nonlinear corrections depend on the incident power 
differently in the two theories. 

In Fig. m we plot the absolute value of the nonlinear 
correction to the particle polarizability normalized by its 
volume as a function of radius for the same values of 
frequency as were used in Fig. |4l We denote the quantity 
being plotted by ^nl and 



as (///at) , in the "quantum" case 
03 \////at , in the "classical" case . 



(53) 



The nonlinear effects should be observable in measure- 
ments with incoherent light if |^nl| ?i 1- If |^nl| ^ 1, 
detection of the nonlinear effects requires coherent laser 
excitation and utilization of the effect of synchronism. 

The parameters used in Fig. |5] are such that the ap- 
proximate analytical formulas for xi p4)) and xa p7|) 
are fairly accurate, as was demonstrated above. Cor- 
respondingly, we have used these formulas to generate 
the curves, which are displayed in Fig. [B] To obtain the 
"quantum" curves, the following procedure has been fol- 
lowed. First, we have computed the function X3(a) ac- 
cording to (j37p for each frequency considered. Then we 
have computed Xii*^) according to for the same fre- 
quencies. In Eq. p4p . we have accounted for the depen- 
dence of the relaxation constant 7 on a according to p6p . 
The computed function Xii'^) '^as used to compute the 
linear field enhancement factor /i(a) according to (P5|). 
Finally, we have used the functions /i(a) and XaiO') to 
compute 03(0) according to (^B)) . The result was mul- 
tiplied by ///at ~ 4.3 • 10"^"'^. In the "classical" case, 
/Ja was computed according to (j52p , where the relaxation 
constant 7 was size-corrected according to ([55)1 . 

We now discuss the curves shown in Fig. [5] in more 
detail. First, in the "quantum" case, .^nl exhibits an 
unlimited growth with a when a — >■ 00. In the classi- 
cal case, this growth is suppressed. As can be seen, the 
"classical" ^^nl decreases with a in the case uj — O.lojp 
and seems to reach a finite limit in the case uj — ujp/^/3. 
In reality, however, the "classical" curve in Fig. [6jb), 
reaches a maximum at a « 44nm and then slowly ap- 
proaches zero (the range of radiuses, which is necessary 
to see this behavior clearly, is not shown in the figure). 
In the classical theory, the nonlinearity is an effect of the 
finite size, which vanishes in the limit a — > 00. 

Second, when uj = O.lwp, the local-field enhancement 
factor in the "quantum" theory is ~ 10~*. That 
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FIG. 6: (color online) Absolute value of the nonlinear correc- 
tion to the nanoparticle polarizability, i^NL, computed using 
Eq. (|50|) (the "quantum" curves) and Eq. (|5ip (the "classi- 
cal" curves) for cj — O.lcjp (a) and cj — Up/y/i ~ 0.58cjp 
(b). Here i^NL = Q3(///at) (the "quantum" curves) and 
i^NL = Psyljlzt (the "classical" curves). The incident power 
is / = lOkW/cm^, J/J^t « 4.3 • 10-". To compute 03, the 
internal field enhancement factor has been taken into account 
according to H46p . 



is, the field is effectively screened in the interior of the 
nanoparticle. Correspondingly, the nonlinear effect is 
very weak. In the "classical" theory, the field enhance- 
ment factor is different, namely, it is |/3ip ^ 1. This 
dramatic difference is explained by the fact that the clas- 
sical theory considers the induced electron density near 
the nanoparticle surface where the electric field is not en- 
tirely screened. At the Frohlich frequency, to — uip/^/S, 
the situation is quite different; we have 



7-10" 



Correspondingly, the "quantum" nonlinear correction 
can become very large at the Frohlich frequency; this is 
illustrated in Fig. |6l[b) . This result is probably unphys- 
ical - one can not expect that ^nl ~ 10"^ at the mod- 
est incident intensity of lOkW/cm^. The classical curve, 
however, is still bounded at the Frohlich frequency be- 
low 10""'. We can conclude therefore that the local- field 
correction plays a disproportionate role in the quantum 
theory and that, if used unscrupulously, it can predict 
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utterly unrealistic magnitudes of the nonlinear effect. 

VI. SUMMARY OF FINDINGS AND 
DISCUSSION 

In this paper, we have further developed the quantum 
theory of Refs. HQ (the HRFR model). The goal was 
to describe the frequency and size dependence of linear 
and nonlinear optical susceptibilities of spherical metal 
nanoparticles. We have used the HRFR model without 
modification but have managed to simplify the previously 
published expressions to a point where these expressions 
became amenable to direct numerical implementation. 
Then, we have computed the linear and nonlinear sus- 
ceptibilities numerically for various frequencies and var- 
ious particle sizes and compared the obtained results to 
Rautian's analytical approximations. Previously, numer- 
ical computations of this kind have been hindered by the 
overwhelming computational complexity of the problem. 
We have also compared the predictions of the quantum 
theory of size-dependent optical susceptibilities with the 
predictions of a purely classical theory of Ref. [l^ The 
following findings can be reported: 

1. We have found that the approximate formulas de- 
rived by Rautianl are surprisingly accurate in a 
wide range of parameters despite the use of a num- 
ber of approximations. In particular, we have, for 
the first time, verified from first principles the cor- 
rectness of the widely-used finite-size correction to 
the Drude relaxation constant (pB)) . 

2. We have found that, for sufficiently small values 
of radius and frequency, Rautian's approximations 
break down due to the discreteness of electron en- 
ergy levels. At sufficiently small frequencies, a sil- 
ver particle with a < 4nm in radius behaves as a 
dielectric. This effect is illustrated in Fig. [l{a) for 
a = 2nm. 

3. We have found that phenomenologically account- 
ing for the local-field correction (see Sec. |IV] for 
details) does not remove the two main difficulties, 
which are encountered in the HRFR model, namely, 
the incorrect small-w asymptote for the linear sus- 
ceptibility xi and the absence of a "bulk" limit for 
the nonlinear susceptibility X3. It appears that ob- 
taining the correct asymptotes requires the rigorous 
account for the Hartree interaction potential. It is 
also conceivable that obtaining the correct large- a 
asymptote requires accounting for the retardation 



effects. However, the classical theory of Ref. [H is 
quasistatic but does not possess a large-a anomaly. 
This suggests that the main focus in further devel- 
opment of Rautian's theory should be on a more 
accurate inclusion of Hartree interaction. 

One additional comment on the theory devel- 
oped here are necessary. First, we have computed 
only a particular case of the nonlinear susceptibility 
X^'^^(w; aji,cj2,W3). More specifically, the coefficient 
X3 defined in ([27)) is related to the latter quantity 
by Xs = ^a^x'^^H^T^^ However, considera- 

tion of transient processes, generation of combination 
frequencies and harmonics requires the knowledge of 
x'''^'(w; wi, W2,W3) as a function of all of its arguments. 
This is an important consideration. High incident inten- 
sities are usually obtained in short laser pulses. More- 
over, many modern photonics applications such as waveg- 
uiding, etc., utilize short wave-packets. Therefore, a 
proper description of optical nonlinearities in a transient 
process is very important. Generalizing the mathemat- 
ical formalism described in this work to include three 
independent frequencies is not conceptually difficult, al- 
though can lead to cumbersome calculations. 

In summary, the HRFR model forms a perfect theo- 
retical framework for studying optical nonlinearities and 
finite-size effects in nanoparticles. The only viable al- 
ternative to using this model is to resort to density- 
functional theory (DFT). In a recent paper—, we have 
applied DFT to study the nonlinear electromagnetic re- 
sponse of metal nanofilms, but only at very low frequen- 
cies, well below plasmonic resonance of the system, and 
neglecting the relaxation phenomena. Higher frequen- 
cies, which are of interest in plasmonics, can be stud- 
ied with the use of time-dependent DFT (TDDFT). Al- 
though TDDFT has been used successfully to compute 
linear response of nanoparticles and, in particu- 
lar, to study the effects of surface adsorption of various 
molecules on the relaxation phenomena in metal^^i^i, the 
difficulties here are formidable. Most importantly, there 
is almost no hope of obtaining analytical approximations 
within DFT. It appears, therefore, that devising a way 
to include the Hartree interaction potential in the master 
equation p3p would be a useful and practically-relevant 
development of the HRFR model and of Rautian's theory. 
Perhaps, some elements of DFT can be used to achieve 
this. 
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Appendix A: Functions (?i(«;) and gain) 



The integrals ([55)1 and ([55)1 can be evaluated analyti- 
cally with the following results: 



2 6k 

-9 + -- 
X (35 -I- 5k 



6k2 
— 6k^ 



2k3 

-f 8k3. 



2VT 



315 
16k4) 



(Al) 



640Kg3(K) = [128 + k{k + 2)(88 + 5k(3k - 8))] VTT^ 



[128 + k(k - 2)(168 - 5k(3k + 8))] y/l - k 

(A2) 



-15k'' In 



1 + 



1 + 



The above expressions have been obtained from p3p and 
psp without using any approximations. However, it 
should be kept in mind that ((33|) is valid for —1 < k, 
while (pS)) is valid for — 1 < k < 1. Since k = hu/Ep, 
we are interested only in the region k > 0. The functions 
(7i(k) and 93 (k) in the interval < k < 1 are shown in 

Fig.m 



